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Abstract 

As showed in (Fiedler, 1990), any polynomial can be expressed as a characteristic poly¬ 
nomial of a complex symmetric arrowhead matrix. This expression is not nnique. If the 
polynomial is real with only real distinct roots, the matrix can be chosen real. By using 
accurate forward stable algorithm for computing eigenvalues of real symmetric arrowhead 
matrices from (Jakovcevic Stor, Slapnicar, Barlow, 2015), we derive a forward stable algo¬ 
rithm for computation of roots of such polynomials in O(n^) operations. The algorithm 
computes each root to almost full accuracy. In some cases, the algorithm invokes extended 
precision routines, but only in the non-iterative part. Our examples include numerically 
difficult problems, like the well-known Wilkinson’s polynomials. Our algorithm compares 
favourably to other method for polynomial root-finding, like MPSolve or Newton’s method. 

Keywords Roots of polynomials; Generalized companion matrix; Eigenvalue decomposition; 
Arrowhead matrix; High relative accuracy; Forward stability 

MSC 65F15, 65G50, 15-04, 15B99 

1 Introduction and Preliminaries 

Polynomials appear in many areas of scientihc computing and engineering. Developing fast al¬ 
gorithms and reliable implementations of polynomial solvers are of challenging interest. Famous 
example by James H. Wilkinson in 1963 m, usually referred to as Wilkinson’s polynomial, is 
often used to illustrate difficulties when finding the roots of a polynomial. The polynomial of 
order n is defined by a simple formula: 

20 

tFn (a;) = — i) = {x — 1) {x — 2) ■ ■ ■ {x — n). 

i=\ 

For example, the location of the roots of W 20 is very sensitive to perturbations in the coeffi¬ 
cients, so that in m, Wilkinson said: ’’Speaking for myself, I regard it as the most traumatic 
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experience in my career as a numerical analyst.” Many methods for finding roots of polynomials 
with ever increasing accuracy have been developed since (see for example m, IS]). 

In a, Miroslav Fiedler showed that any polynomial can be expressed as a characteristic 
polynomial of a complex symmetric arrowhead matrix. This expression is not unique. If the 
polynomial is real with only real distinct roots, the matrix can be chosen real. We have the 
following theorem: 


Theorem 1. Theorem 3] Let u (x) be a polynomial of degree n, 


u (x) = x" + px" 

"^ + r (x), 

(1) 

Let 

D = diag((ii,.. 

where dj are all distinct and u{dj) / 0. Let 

• ; ^n—1)7 

(2) 


n—1 

(a;) =n ’ 

n—1 

a = -p-'^dj, (3) 

^=[Cl C2 ••• Cn-l]^, 


where 


—u{dj) 

v'{dj) 


Then the symmetric arrowhead matrix 


—u{dj) 

n—1 

\[{d,-d,) 

i=l 


(4) 


A = 


D z 

T 

z a 


has characteristic polynomial {—!)"'u{x). 


(5) 


If u (x) has only real distinct roots and the dj’s interlace them, then A is real. 

Fiedler concludes his paper by stating ’’One can hope to obtain, by some sophisticated special 
choice of the numbers dj, stable or even universal algorithms for solving algebraic equations. ”0 
In [8] the authors developed a forward stable algorithm for computing eigendecomposition of 
a real symmetric irreducible arrowhead matrix, which is exactly the matrix A given by Theorem 

^In a report by Corless and Litt [2], the matrix A from theorem [T] is referred to as generalized companion 
matrix not expressed in monomial basis. In this case, the basis is the Newton basis. 
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!□ More precisely, the algorithm from [5] computes each eigenvalue and all individual compo¬ 
nents of the corresponding eigenvector of a given arrowhead matrix of floating-point numbers 
to almost full accuracy in 0(n) floating point-operations, a feature which no other method has. 

In this case, we are interested only in the roots of u, that is, in the eigenvalues of A from (jSjl, 
each of which is computed independently of the others in 0{n) operations. This, together with 
independent computation of elements of z, makes our algorithm suitable for parallel computing. 

In this paper, we propose a new two-step algorithm: given a polynomial u of the form ([T]) 
whose coefficients are given floating-point numbers, 

1. compute the generalized companion matrix A from ([5]), where the elements of z and a 
need to be computed in double the working precision, and then 

2. compute the roots of u as the eigenvalues of A by using modified version of the forward 
stable algorithm aheig from [8l Algorithm 5]. 

The organization of the paper is the following. In Section [2l we describe our algorithm 
named poly-aheig (POLYnomial roots via ArrowHead EIGenvalues). In Section [3l we analyse 
the accuracy of the algorithm and give forward error bounds - in Section 13.11 we analyse the 
accuracy of the computed matrix A, and in Section [321 we analyse the accuracy of the computed 
inverse of the shifted matrix A. In Section 13.31 we discuss possible ways to find the diagonal 
elements of the matrix D which interpolate the roots of u. In Section 13.41 we discuss different 
implementations of the double the working precision, including extended precision routines from 
[3] and the Compensated Horner scheme from [H Algorithm 4]. Finally, in Section 01 we illustrate 
our algorithm with two numerically demanding examples and compare it to the methods from 
[I] and [5]. 


2 The algorithm 

The eigenvalues of the arrowhead matrix A from ([5l) are the zeros of the function 

PA (A) = a — \ — {D — \I)~^ z. 

The forward stable algorithm for solving EVP of arrowhead matrices [8] computes all eigen¬ 
values to almost full accuracy. The algorithm is based on shift-and-invert strategy. Let di be 
the pole which is nearest to A. Let Ai be the shifted matrix. 


Aj = A — dji = 


Di 

0 

0 

Zi 

0 

0 

0 

Ci 

0 

0 

T>2 

Z2 

T 

^1 

Ci 

T 

^2 

a 


( 6 ) 


^In [8], the arrowhead matrix A is called “irreducible” if dj are all distinct and Zj ^ 0, j — 1,... ,n — 1. 
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where 


Di = diag(di - di,..., di-i - di), 

D 2 = diag(cii+i - di,..., dn-i - di), 

= [ Cl C 2 • • • C*-i (7) 

Z2 = [ Ci+1 Ci+2 • • • Cn-1 ]^, 

a = a — di. 


Then, 


A — — di, 

V 


where v is either largest or smallest (first or last) eigenvalue of the matrix 

-1-1 


where 


A-^^{A-da)-^== 


DI 

wf 

0 

0 


Wl 

b 


0 


W 2 
W2 

0 


1/0 


0 

1/0 

0 

0 


Wl = -D^ 

W2 = -D2^Z2y, 
Si 


1 


b = 72 (-“ + zi zi + Z 2 D 2 Z 2 ). 


( 8 ) 


(9) 


Notice that all elements of the matrix A~^ are computed with high relative accuracy, except 
that in some cases the element b needs to be computed in double the working precision (for 
details see m)- Also, the elements of z (the Horner scheme) and a (the trace preservation 
formula) of A need to be computed in double the working precision. Notice that our algorithm 
requires computation in higher precision only in the finite part, unlike algorithms from HIS], 
which require usage of higher precision in the iterative part. 

The described procedure is implemented in the algorithm poly-aheig. 
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Algorithm 1 
A = poly_aheig(?/, D) 

% Computes the roots A of the polynomial u(x) from ([1]) of order n, with 
% distinct real roots. D is defined by ([2]) and its entries interlace the roots 
% of u{x), see Section [331 for details. 

% Compute the values of u{x) in the interpolating points dj using double 
% the working precision, 
for j = 1 : n — 1 

^doubleiJ) ^(^(j)) 

end 

% Compute vector v from Theorem [1] using double the working precision, 
for j = 1 : n — 1 

VdoubieU) = - d{l:j-l,j + l:n- 1)) 

end 

% compute a from Theorem [1] using double the working precision. 

n—\ 

(^double P dj 

% compute vector z from Theorem [1] using double the working precision, 
for j = 1 : n — 1 

CdoubleiJ) ~ 'J Sdoubleij)/Vdouble (j) 

end 

% call modified algorithm aheig 
for k = 1 : n 

X{k) = aheig_mod(T>, Zdoubie, adoubie, k) 

end 
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Remark 1. The algorithm aheigjmod is a simple modification of the algorithm aheig from 
[U Algorithm 5]. The accuracy of the algorithm aheig is essentially based on the assumption 
that all elements of the matrix A~^ from ([5]) can be computed with high relative accuracy, 
that is, = [A~^]ji{l + KjiSM), for some modest Kji. For all elements of A~^ but b, 

this accuracy is achieved by computing them in standard precision using the standard precision 
copies of z and a. If, according to the theory from [8], the element b needs to be evaluated in 
double the working precision, formula ([9]) is evaluated using Zdoubie and Odoubie in order to obtain 
full possible accuracy. The details of the analysis follow. 


3 Accuracy of the algorithm 


The error analysis of the algorithm aheig is given in [8l Sections 3 and 4]. This analysis assumes 
that A is the given matrix of floating-point numbers. Here, however, A is computed by using 
formulas (HED, which must be taken into account. We assume that computations are performed 
either in the standard floating-point arithmetic with the machine precision em = ~ 1.1102- 

10“^® (see P Chapter 2] for details) or with double the working precision with the machine 
precision « 1.2326 • 10“^^. H 

Let us first consider the errors in the polynomial evaluation. The standard method for 
evaluating polynomial u{x) is the Horner’s method P Section 5.1]. Let 

n 

u(x) = ^ ao = 1, (10) 

i=0 


and let 


cond{u,x) = 


n 

Wi\ kl 

i=0 


E 

i=0 


u{x) 

u{x) 


( 11 ) 


Notice that cond{u,x) > 1. Let Horner{x,u) denote the value of u(x) computed in floating 
point accuracy by the Horner scheme. Then, the relative error in the computed value is bounded 
by P Section 5.1] 0 

\u(x) — Horner (u,x)\ ,, , 

-j—-< cond(u,x) X znSM- 

\u (x)| 

Thus, when Horner{u,x) is evaluated in double the working precision, the relative error is 
bounded by 


\u{x) - Hornerdoubieiu,x)\ 


< cond{u,x) X 2nej^. 


Therefore, 


\u{x)\ 

Horner doubie{u,x) = {I + Kxe\j)u{x), 


( 12 ) 


®Thus, the floating-point numbers have approximately 16 significant decimal digits. The term “double the 
working precision” means that the computations are performed with numbers having approximately 32 significant 
decimal digits, or with the machine precision equal to e%[. 

'^In [ni[S], the bounds are expressed in terms of quantities 7 fe = For the sake of simplicity, we use 

standard first order approximations yt ~ ksM- 
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where 


( 13 ) 


< cond{u,x) X 2n. 

Notice that, if cond{u, x) is uniformly bounded, 

cond{u,x) < —, (14) 

£m 

then 

\nx\ < 2n. (15) 

Other two options to obtain bounds similar to (jl2ll3p are to evaluate all parts of the respec¬ 
tive formulas by using extended precision routines from [3] or compensated Horner scheme from 
O Algorithm 4] (see Section [33] for details). 

We now consider the accuracy of the computed matrices A, Aj and A~^ from ([5]), ([ 6 |) and 

©• 


3.1 Accuracy of A 

Let A denote the matrix A computed according to Algorithm 1 , 


A = 


D zW 


Here and are computed in double the working precision which we denote by superscript 
{d). Let 




Ad) Ad) 


Ad) 
i’n—1 


1 T 


By combining ([4]) and (I12p . the standard first order error analysis in double the working 
precision, gives 




-u{dj){l + 


n—1 


■(1 + ^3)(l + £ 4 ), 


(16) 


{dj - di){l -h ei)(l -b (n - 3) 62 ) 

A i=l 


where |ei, 2 , 3 , 4 | < Therefore, 


^ “ 0(1 + 


(17) 


where, by using (ITsp . 


Ad) 


I T (^ 1) u-|-l 

< ---h 1 < n • cond{u, dj) --—. 


Similarly, applying the standard first order error analysis in double the working precision to (|3]) , 
gives 

aA) = Q,(^i _|_ 
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where 




< 


n—l 

l®i| + Z] l^il 
_ 

\a\ 


-(n - 1) = Ka{n - 1). 


(18) 


3.2 Accuracy of ^ 

Let A~^ denote the matrix A~^ computed according to Algorithm 1 from the matrix A. All 
elements of A~^ but possibly b, are computed in standard precision using the standard precision 

copies of and Let Q and a denote ^ and rounded to the nearest standard 

precision number, respectively. Then 


Cj = 0 (l + kq£m) , j = 

d = cr (1 T ) 


where, by using (fT7)l ~ (fT8]) . 


kc.l ^ ( 


i^di I + (?^ - 1) 


Em + 1, j = l,...,n-l. 


l^al < Ka{n - 1 )£m + 1- 


Further, according to (fT3l) - (fT^ . if 

cond{u,dj) <—, j = l,...,n —1, 
£m 


(19) 

( 20 ) 


( 21 ) 


then (fTOl) holds with 

<n + 2, j = 1,.. 

. ,n - 1, 

(22) 

and if 

Ka < 


(23) 




then (1201) holds with 

\Ka\ <n. 


(24) 


For j ^ {i,n}, similarly as in [8l Proof of Theorem 4], the standard first order error analysis 
gives 

■'""I - 

Similarly, assuming that ([211) and ([2^ hold, for j ^ {i, n} we have 


= fi 

“0 


“0(1 + 


{dj di)(^i 


{dj — di)Q{l + kqEm) 

(1 + \^ji\ — (2^ T I”)- 










Finally, 


[A'l- = fm7%) = n 


Ci(l + 


1 


— ^ (1 “1“ |^m| ^ “t“ 3). 

si 

We now analyze the accuracy of the computed element b. Let 

|q:| + \di\ + \z^D^ ^Zi\ + \z2 D 2 ^"^ 2 ! 


Kh = 


— a + zfD-^ ^zi + Z 2 D 2 ^Z 2 \ 


(25) 


where Di, D 2 , zi, Z 2 and a are defined by ([7]). We have two cases. First, if 


1 , 

then b is computed in standard precision using Q and d. Let b denote the computed b. The 
standard first order error analysis of ([9]) gives 

a(l + KaEM) — di + — - 

i=i 
j¥=i 

= 6(1 + KbEM), 


b = fl 


Q (1 + 


+ K-QEmY 

dj — di 


where 


\i^b\ < (n + 2 + max{2max |ka-|, |kq|}) • Kb + 2|ka| + 3. 


Additionally, if (l2T]l and (1^ hold, then (1^ and (IMl) hold, as well, and 

|Kft| < (3n + 6) • iLfe + 2n + 7. 


Second, if 

Kb > 1, 

then, according to the theory from [8], the element 6 needs to be computed in double the working 

precision using (^j and in order to obtain full possible accuracy. The standard hrst order 
error analysis of Q in double the working precision gives 


= fi 


Q(1 + 


n —1 a2('i 

i=i 


dj — di 


— 6(1 + 


where 


Mi 


< (n + 2 + max{2 max I}) • Kb + 2\kY’ I + 3 
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Finally, let 


= max{2max (26) 

Ai jz^i 

If, in addition to (f2Tl) and ([23]), 

Kb <—, (27) 

£m 

and 

• iffe < —, (28) 

Em 

then 

b = = 5(1 + Kbeiii), 

where 

l^fel < n + 5. 

The above results are summarized in the following lemma: 

Lemma 1. Let fill) and (E^ hold, and let Kb he defined by l[25\) . For all non-zero elements 
of the matrix A~^ from |3) computed according to Algorithm 1 and Remark 1, except for the 
element {A~^]ii, we have 

Ai ^\kl = Ai ^]fcz(l + l^kFAl), \l'^kl\ ^ (2u + 7). 

For the computed element b = we have the following: if Kb ^ 1, then 

5 = 5(1 + KbEM), l^^bl < (3n + 6) • iFfe + 2n + 7. 

// iFfe 3> 1 and if (2l\) and hold, then 

5 = 5(1 + Kfee^), |Kb|<n + 5. 


The forward error of the computed roots is bounded as follows: 
Theorem 2. Let [2l\) and (2^ hold, and let Kb he defined by E5l) . Let 

A = A(1 + k\£m) 


be the root of u{x) computed according to Algorithm 1 and Remark 1. If Kb 1, then 
|ka| < 3-v/n[(3n + 6) • iF;, + 2n + 7] + 3.18n (yy/n + l) + 4, 
and if Kb I and and hold, then 

I^^aI < (6n + 21)AAj' + 3.18n (\/n + l) + 4. 


Proof. Using the same notation as in [8l §3], the first summand in the above bound for k\ follows 
from [HI Theorems 5 and 6], while the second summand is the error bound for bisection from 
[ini §3.1], □ 
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3.3 Choosing dj 

Finding values of dj which interpolate roots is not an easy task. Articles dealing with computing 
roots of polynomials usually assume that the initial approximations of the roots are known (see 
0 , m)- Another approach, used in [T], is to define the polynomial neighborhood of u{x) as 
the set of all polynomials with coefficients having din common digits with the corresponding 
coefficients of u{x), where din is predefined input precision. Then, the root neighborhood is the 
set of the roots of all polynomials in the polynomial neighborhood of u{x). 

Since our polynomial is real with only real distinct root our proposal is simpler. Here are 
some heuristics: 

let u{x) be the reverse polynomial of the polynomial u{x) from cni), 

u{x) = x^u{l/x) = qqx^ + aix'^~^ + + • • • + a „_ 2X ^ + an-ix + 1 . 

Since the roots of u{x) are the reciprocals of the roots of u{x), we have two options for the 
values dj\ 

- use the roots of u'{x), or 

- use the reciprocals of the roots of u'{x). 

Depending on the magnitude of the roots, their distribution and relative gaps, one of the 
methods, or a combination, is expected to work, see Section 0] for examples. 

3.4 Implementation of double the working precision 

We tested three different implementations of double the working precision: 

- convert all quantities to variable precision by Matlab command sym with parameter ’ f ’, 
and then evaluate the respective formulas - this is 300 to 1000 times slower than standard 
precision. 

- convert all quantities from standard 64 bit REAL(8) to 128 bit REAL(16) in Intel ifort 
[7], and then evaluate the respective formulas - this is only 3 times slower, 

- evaluate respective formulas by using extended precision routines add2, sub2, mul2, div2, 
and sqrt2 from [3] - this is 0(10) times slower. In these routines double the working 
precision is simulated by keeping each number as a pair consisting of higher and lower 
part of mantissa. For example, let 

[z, zz] = add2{x, xx, y, yy) 

where all quantities are floating-point numbers with t binary-digits mantissa. Then 

jz + zz - [(x -t- xx) -h (y-h yy)]\ < {\x + xx\ + \y + yy\)2~^^^~^'>. 

If XX = 0 and yy = 0, then (exactly) z + zz = x + y. We see that this is nearly equivalent 
to using double the working precision (the precision is instead of e^). 

The evaluation of the polynomial u{x) can also be successfully performed by Compensated 
Horner scheme from [5l Algorithm 4], where both quantities h and c from this algorithm must 
be preserved for subsequent computations by extended precision routines. 
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4 Numerical Examples 


We illustrate our algorithm with two numerically demanding examples. Here double the working 
precision in Algorithm 1 was implemented with extended precision routines from [3]. 

Example 1. The coefficients of Wilkinson polynomial Wis are, row-wise 


1 

-662796 

10246937272 

-14710753408923 

2353125040549984 

-30321254007719424 

6402373705728000 


-171 

22323822 

-147560703732 

102417740732658 

-7551527592063024 

34012249593822720 


13566 

-549789282 

1661573386473 

-557921681547048 

17950712280921504 

-22376988058521600 


In this example the interpolating points dj can be computed by both ways described in 
Section (3^ as roots of u'{x) or as the reciprocals of the roots of u'{x). For example, in the 
latter case we have 


maxATft = 214.5 1, max{cond(n, dj)} = 2.62 • 10^^, = 26.8, 

j 

so by Theorem m the roots of Wis are computed by Algorithm 1 to (almost) full accuracy, in a 
forward stable manner. 

The roots computed by Matlab [9] routine roots, MPSolve [T] (with 16 decimal digits), 
Algorithm [1] and Mathematica |14] with 100 digits of precision (properly rounded to 16 decimal 
digits), are, respectively: 


^{roots) 

^{MPSolve) 

^{poly^aheig,Math) 

18.00001193040660 

18.00000000000000 

18 

16.99987506992020 

16.99999999999993 

17 

16.00057853967064 

15.99999999999455 

16 

14.99841877954789 

15.00000000000043 

15 

14.00282666587300 

13.99999999999777 

14 

12.99649084561071 

12.99999999999819 

13 

12.00308090986650 

12.00000000000329 

12 

10.99809154207482 

11.00000000000163 

11 

10.00081885564820 

9.999999999998594 

10 

8.999776556759201 

9.000000000000055 

9 

8.000029075840132 

7.999999999999923 

8 

7.000002735870642 

7.000000000000000 

7 

5.999998227088450 

5.999999999999999 

6 

5.000000283698958 

5.000000000000000 

5 

3.999999981972712 

4.000000000000000 

4 

3.000000000132610 

3.000000000000000 

3 

2.000000000018936 

2.000000000000000 

2 

0.999999999999808 

1.000000000000000 

1 


®We use Wi8 since all its coefficients are exactly stored as 64-bit floating-point numbers. 
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j 

dj 

cond{u, dj) 

1 

5.277655813324802e+13 

4 

2 

1.759218604441599e+13 

3.58 • 10^6 

3 

6.253878705847983e-16 

12.4 

4 

2.627905491153268e-16 

46.4 


Table 1: Interpolating points dj and cond{u,dj). 


Since for every root, the corresponding quantity ^ 1, the algorithm polyjaheig computes 
fully accurate roots, using only standard working precision to compute the corresponding matrix 
A~^ and its absolutely largest eigenvalue. 

MPSolve requires input to be defined as integers. Also, MPSolve uses 21 decimal digits to 
guarantee and obtain relative accuracy of 10“^^, and it uses 234 decimal digits to guarantee and 
obtain 30 accurate digits. 

The Accurate Newton’s method from [5l Algorithm 6] also computes the roots of ITig to 
full accuracy as described in [U Theorem 6]. However, the starting points xq which satisfy the 
conditions of [5l Theorem 6], must be chosen with greater care and must be relatively close to 
the desired root (for example, xq = 17.1 to obtain A 2 = 17, or xq = 1.1 to obtain Aig = 1. Since 
the Accurate Newton’s method takes on average 6 steps to convergence for each root, it needs 
approximately 12n^ effective extended precision computations, while our algorithm needs in this 
case 5n^ extended precision computations to compute the matrix A. 

The results for IT 20 are similar. 

Example 2. Consider the polynomial u of degree 5 with the coefficients 

l.OOOOOOOOOOOOOOOe + 00 
-2.028240960365167e + 31 
7.136238463529799e +44 
-6.277101735386680e + 57 
4.1813897247244916 + 42 
-6.1897001964269006 + 26 


or symCu,’f’) 


1 

-20282409603651670423947251286016 

713623846352979940529142984724747568191373312 

-6277101735386680066937501969125693243111159424202737451008 

4181389724724490601097907890741292883247104 

-618970019642690000010608640 

In this example the interpolating points dj are efficiently computed as the reciprocals of the 
roots of u'(x). The values dj and cond(u, dj) from ifTTI) are given in Tabled) 

For the decreasingly ordered roots of u, Xk, k = 1,2, 3,4, 5, the corresponding quantities Kb 
from ([ 25 ]), from (l26)) and their respective products from (l28]l . all rounded up, are given in 
Table [2l 


13 





k 

Kb 

Ad) 


1 

1 

3.6 • 10^^ 

3.6 • 10^^ 

2 

3.01 • 10^® 

4.7- 10^ 

1.42 • 10^® 

3 

3.01 • 10^5 

4.7- 102 

1.42 • 10^® 

4 

12.6 

3.58 • 10^^ 

4.48 • 10^® 

5 

12.6 

3.58 • 10^^ 

4.48 • 10i8_ 


Table 2: Values Kb, and ■ Kb- 


We see that the condition ([27]) is always fulfilled. Also, = 1 from (fTB|l . so (f25|l is fulfilled. 
The condition (|28p does not hold literally. However, we have • Kb ^ which is sufficient 
to obtain almost full accuracy. 

The roots computed by Matlab [9] routine roots, MPSolve [T] (with 16 decimal digits). 
Algorithm [1] and Mathematica [H] with 100 digits of precision (properly rounded to 16 decimal 
digits), are, respectively: 


^{roots) 


^{MPSolve) 


^{poly^aheig,Math) 


2.028240960365167e + 31 2.028240960365167e + 31 

1.759218604441600e + 13 + 1.538e + 8i 1.759218604441608e + 13 
1.759218604441600e + 13 - 1.538e + 8i 1.759218604441591e + 13 

0 4.440892098500624e-16 

0 2.220446049250314e-16 


2.028240960365167e+31 
1.759218623050247e+ 13 
1.759218585832953e+ 13 
4.440892098500624e- 16 
2.220446049250314e- 16 


We see that the roots computed by Algorithm [T] coincide fully with roots computed by 
Mathematica. Here, in addition to the elements and of the matrix A, the element b of 
^2 ^ was computed in double the working precision. 

Again, MPSolve requires input to be defined as integers, and it uses 21 decimal digits to 
guarantee and obtain relative accuracy of 10“^^, and uses 234 decimal digits to guarantee 30 
accurate digits. 

Here the Accurate Newton’s method from [5l Algorithm 6] also computes the roots to full 
accuracy, provided the respective starting points are chosen with greater care. However, the 
conditions of O Theorem 6] cannot be used - for example, for the largest root Ai, there is no 
starting point xq which satisfies the conditions, except Ai itself. For A 2 , the starting point xq 
which satisfies the conditions can differ from A 2 in just last digit. 
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